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ABSTRACT 


The  use  of  direct  counts  of  molecular  energy  level  degeneracies,  P(E), 
at  internal  energy,  E,  in  RRKM  unimolecular  reaction  rate  calculations  is 
clarified.  Highly  accurate  densities  of  states  for  use  in  practical  rate 
calculations  can  be  obtained  by  utilizing  an  energy  level  "grain"  size  of 
1-5  cm-1  and  by  suitable  averaging  of  P(E)  over  small  energy  intervals.  An 
easily  programmed  algorithm  for  rounding  energy  levels  to  multiples  of  this 
"grain"  size  is  given  and  applied  to  sum  and  density  calculations. 


1. 


INTRODUCTION 


The  advent  of  highly  efficient  algorithms  [1-3]  for  direct  computation 
of  energy  level  degeneracies,  P(E),  and  sums,  ZP(E),  now  permit  such  computa- 
tions to  be  performed  with  great  accuracy  over  the  entire  energy  range  of 
chemical  interest.  In  practical  calculations,  these  algorithms  typically 
require  5000-20,000  computer  storage  spaces  and  an  execution  time  of 
0.1-100  sec. 

Using  energy  levels  or  fundamental  harmonic  oscillator  frequencies  of 
individual  degrees  of  freedom  as  input,  direct  counting  algorithms  yield 
a linear  array,  P't  whose  elements  P'(I)  contain  the  number  of  different  ways 
of  arranging  the  energy  (1-1)  G in  the  available  energy  levels,  where  G is 
a chosen  "grain"  size.  For  present  purposes,  P(E)  = P’(R(E/G)  +1),  where 

the  quantity  R(E/G)  denotes  the  integer  nearest  to  the  value  (E/G).  Thus, 
P(E)  is  effectively  the  number  of  levels  around  energy  E within  an  energy 
interval  of  length  G.  The  only  errors  in  P(E)  that  may  arise  using  this 
method  result  from  the  necessity  of  altering  the  input  energy  levels  (or 
fundamental  frequencies)  to  become  integral  multiples  of  G.  An  algorithm 
proposed  by  Beyer  and  Swinehart  [1]  (BS),  which  has  been  compared  to  other 
schemes  by  Stein  and  Rabinovitch  [2],  computes  P(E)  and  £P(E)  for  molecular 
degrees  of  freedom  whose  energy  levels  are  evenly  spaced  (harmonic 
oscillators).  For  maximum  accuracy,  G is  typically  chosen  to  be  between  1 
cm-1  and  5 cm-1  which  is  generally  within  the  known  accuracy  of  fundamental 
vibrational  frequencies. 

A related  algorithm  by  Stein  and  Rabinovitch  (SR)  computes  P(E)  and 
ZP(E)  for  any  degrees  of  freedom  whose  energy  levels  are  known  [3]  (examples 
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are  Morse  anharmonic  oscillators  and  symmetric  hindered  rotors).  The  SR 
algorithm  computes  exact  P(E)  values  for  sets  of  energy  levels,  each  level  of 
which  has  been  rounded  to  an  integral  multiple  of  G.  Errors  in  calculated 
P(E)  are  negligible  (<  1%  error)  for  molecules  of  chemical  interest  with 
choice  of  a grain  of  1-5  cm-1  [3].  For  a computation  involving  both 
harmonic  and  non-harmonic  degrees  of  freedom,  it  is  straightforward  to  apply 
the  BS  algorithm  for  the  harmonic  degrees  of  freedom  and  then  apply  the  SR 
algorithm  or  other  [4]  for  the  non-harmonic  degrees  of  freedom. 

Recently,  Hayward  and  Henry  (HH)  [5]  have  questioned  the  stability  of 
counting  methods  for  accurately  determining  sums  and  densities  of  states  as 
required  for  RRKM  unimolecular  rate  calculations.  They  have  alleged  that 
such  methods  may  be  highly  inaccurate  partly  because  of  the  inherent  non- 
linearity of  sums  and  densities  as  a function  of  energy,  but  primarily  on 
the  grounds  that  that  lack  of  smoothness  of  P(E)  versus  E can  yield  erroneous 
values  for  densities,  N(E).  It  will  be  shown  that  such  objections  are  not 
valid  for  calculations  using  a small  enough  step  size  with  realistic  vibrational 
frequencies;  for  practical  computations  the  use  of  P(E),  derived  from  direct 
counting  methods,  to  obtain  N(E)  is  highly  accurate. 

In  this  paper  we  shall  discuss  the  correct  application  of  the  output  of 
these  computations,  P(E),  for  evaluating  densities  and  sums  of  states. 

Attention  is  also  given  to  the  method  used  for  rounding  energy  levels  and 
frequencies  to  multiples  of  G. 

2.  GENERAL  DISCUSSION 

In  order  to  obtain  an  overall  rate  constant  for  an  elementary  reaction 
occuring  over  a range  of  energies,  rate  constant  elements  are  individually 
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evaluated  over  small  energy  intervals,  from  E to  E + £,  and  then  summed. 

In  particular,  the  RRKM  unimolecular  rate  constant  for  thermally  induced  reac- 
tions, kuni,  may  be  numerically  evaluated  as  [6] 

n£ 


kUni  - a I 
n=o 


'(  ZP(E+))exp(-n£/kT) 

E+=0 

n£ 

1 + b(  ZPCE^JJ/NCEo  + n£)p 
E+=0 


(1) 


where  a and  b are  constants  which  depend  on  the  specific  reaction  type,  p 

is  the  pressure,  E0  is  the  threshold  energy  for  reaction  and  E*  is  the  active 

energy  in  the  activated  complex.  For  computations  using  direct  counting,  the 
n£ 

term  e_^P(E)  is  obtained  by  summing  all  P'  (I)  elements  up  to  I = (E/G)  + 1. 
The  term  N(E0  + n£)  is  evaluated  by 


N(E0  + n£) 


1 R(4> 

j Z (P'(R  (E0  + n£)/G)  + 1 + M) 

% n 


(2) 


This  expression  simply  sums  all  P(E)  elements  in  the  energy  range 
E0  + (n  - y)£  to  E0  + (n  + -j-)£.  This  formula  (2)  is  strictly  correct  only 
when  £/G  is  not  an  integer  in  which  case  care  must  be  taken  to  count  each 
P'  element  only  once. 

For  the  sum,  given  by  formula  (1),  to  accurately  approximate  an  integral 
quantity,  the  value  £ must  be  chosen  small  enough  to  ensure  a linear  variation 
of  the  integral  quantity  over  the  range  of  each  summation.  Thus,  it  is 
necessary  for  P(E),  IP(E)  and  exp(-E/kT)  to  be  linear  over  each  energy  interval 
of  length  £.  A choice  of  £/cm-1  between  50  and  200  is  small  enough  to  ensure 
such  behavior,  even  for  large  molecules.  To  illustrate,  P(E)  values  derived 
for  a 120  frequency  harmonic  oscillator  are  shown  in  Figure  1 near  energies 
typical  of  unimolecular  chemical  reaction  (10,000  cm-*  and  20,000  cm-1).  For 


any  calculation  in  which  non-linearity  of  this  type  is  suspected,  it  is 
straightforward  to  perform  a second  calculation  with  a reduced  value  of  Z. 
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To  illustrate  the  roughness  of  P(E)  versus  E and  the  supposed  dependence 
of  N(E)  on  G due  to  this  roughness,  HH  chose  a harmonic  oscillator  model  in 
which  all  frequencies,  Vj,  were  multiples  of  100  cm-1.  They  reported  that 
N(E)  is  dependent  on.  G for  G/cm-1  = 10,  50  and  100.  They  have  chosen  a 
pathological  case  which  can  only  result  by  pre-gaining  each  Vj  to  be  exactly 
divisible  by  the  value  nG,  where  n=2,  3,4,  ...  In  such  cases,  (n-1)  out 
of  n members  of  the  P’  array  will  necessarily  be  equal  to  zero.  For  example, 
in  any  model,  if  the  are  multiples  of  100  cm-1,  there  is  no  way  to 
completely  distribute  an  energy  which  is  an  odd  multiple  of  50  cm-1.  Therefore, 
for  the  G/cm-1  = 50  model,  HH  reported  densities  a factor  of  two  greater  than 

for  the  G/cm-1  = 100  and  for  the  G/cm-1  = 10  model  densities  were  reported  to 

be  higher  by  a factor  of  ten. 

Such  a situation  will  not  arise  if  nZ,  n = 2,  3,  . . .,  does  not  divide 

exactly  into  all  v^.  This  rule,  in  practice,  will  always  be  satisfied  when 

using  frequencies  that  have  not  been  arbitrarily  pre-grained  and  Z/G  > 1. 

To  illustrate  this  point,  shown  in  Table  1,  we  have  computed  P(E)  for  G/cm  1 = 
60  and  50,  with  harmonic  oscillator  model  of  HH.  Note  that  the  variation  in 
P(E)  for  the  former  case  is  very  smooth,  whereas  for  G/cm-1  =50  (n  = 2), 
alternate  P(E)  quantities  are  equal  to  zero. 

Use  of  the  SR  algorithm,  which  rounds-off  each  energy  level , is  even  less 
prone  to  large  fluctuations  in  P'  than  the  BS  algorithm. 

Another  point  of  concern  in  using  these  algorithms  for  N(E)  or  EP(E) 
evaluation  is  the  method  used  for  rounding  fundamental  frequencies  or  energy 
levels.  For  example,  if  all  levels  are  rounded  to  a lower  value,  sums  and 
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densities  will  be  somewhat  overcounted.  More  careful  rounding  can  substantially 
reduce  such  errors.  Of  course,  if  a small  enough  G is  used,  the  absolute  error 
will  be  negligible  no  matter  what  rounding  method  is  used.  A useful  rounding 
scheme  is  as  follows:  Consider  a set  of  im  energy  levels  (or  fundamental 

frequencies)  E^,  i = 1,  2,  . . .,  im,  a linear  array  (with  H0  = 0)  and  a 
grain  size  G.  Let  L(X)  denote  the  greatest  integer  less  than  or  equal  to  X and 
let  represent  the  energy  to  which  the  i-th  energy  level  (or  fundamental 

frequency)  is  rounded.  Then,  for  i = 1 to  im: 

H±  = + Ei/G  - L(Ei/G); 

if  % < Eir  = GL(Ei/G) 
if  % > -^,  Eir  = G(L(EA/G)  + 1)  and. 

Hi  - H±  - 1. 

This  scheme  makes  rounding  decisions  for  each  level  based  on  rounding  of 
previous  levels. 

To  illustrate  the  effectiveness  of  this  rounding  scheme,  in  Table  2 is 
given  sums  and  densities  of  states  have  been  computed  using  the  harmonic 
oscillator  model  of  HH.  Exact  results  are  given  in  column  A (C/cm-1  = 100). 

In  model  B,  which  uses  G/cm-1  = 200,  all  frequencies  were  rounded  downward  or 
unchanged.  This  appears  to  be  the  method  used  by  HH.  Sums  given  by  HH  at 
22,400  cm-1  differ  slightly  from  those  in  Table  2 because  the  former  results 
were  given  as  sums  up  to,  but  not  including  24,400  cm-1.  In  this  present 
paper  sums  are  evaluated  up  to  and  including  the  stated  energies,  using  linear 
interpolation,  if  necessary.  Model  C utilizes  the  rounding  scheme  given 
earlier  in  this  paper  and  is  somewhat  more  accurate  than  the  scheme  of  model 
B.  This  model  is  used  here  for  illustrative  purposes  only,  since  for  actual 
rate  calculations  G will  be  much  smaller  than  the  values  used  here,  and 
resulting  errors  will  always  be  very  considerably  smaller  (<  1%  error). 
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Note,  also,  that  errors  in  EP(E)  follow  errors  in  N(E)  very  closely, 
confirming  the  smooth  behavior  of  the  latter  function. 

In  Table  3 are  given  EP(E)  values  using  G/cm-1  = 10  for  three  anharmonic 
models  given  by  HH,  along  with  deviations  in  EP(E)  resulting  from  the  use  of 
G/cm-1  = 20,  30,  and  100.  Here,  again,  these  are  larger  G values  than  used  in 
practice  and  are  employed  for  illustrative  purposes  only.  The  rounding  scheme 
given  earlier  is  used  in  these  computations.  Note  that  errors  in  EP(E)  are 
uniformly  small  for  even  the  larger  G values.  In  addition,  we  find  that  N(E) 
values,  not  given  here,  are  smooth  functions  of  energy  and  changes  of  N(E)  with 
varying  G follow  changes  of  EP(E)  very  closely  as  for  the  harmonic  oscillator 
model  (Table  1).  Spurious  behavior  of  EP(E)  and  N(E)  as  a function  of  energy, 
as  reported  by  HH,  is  not  apparent;  the  peculiar  results  that  they  obtain  must 
be  due  to  errors  in  their  computational  technique  and/or  their  faulty  round tng 
scheme. 

In  summary,  direct  counting  algorithms  can  be  used  to  supply  accu^r'te 
densities  and  sums  for  RRKM  rate  calculations  as  long  as  grain  sizes  are  small 
(1-5  cm-1),  integration  intervals  are  G and  < 200  cm-1,  and  realistic 
vibrational  parameters  are  employed. 
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TABLE  1 


P(E)  for  Harmonic  Oscillator  Model*1 


G/cia  1 = 60  (present  paper) 


F./ cm  1 

P(E)/10 

9960 

0.2087 

10020 

0.2151 

10080 

0.2243 

10140 

0.2363 

10200 

0.2499 

19980 

37.38 

20040 

38.32 

20100 

39.28 

20220 

41.29 

20280 

42.30 

G/cnf 1 = 50  (HH) 


E/cm  1 

P(E)/10 

9950 

0.0 

10000 

0.3720 

10050 

0.0 

10100 

0.3694 

10150 

0.0 

10200 

0.4415 

19950 

0.0 

20000 

66.03 

20050 

0.0 

20100 

68.05 

20150 

0.0 

20200 

72.96 

aFundamental  f requencies/cm  1 (degeneracies):  3000(2),  1500(2),  1000(2), 

900(2),  500(2),  400(1). 


bFrequencies/cm- 1 after  graining:  3000(2),  1500(3),  1020(1),  960(1), 

900(2),  540(1),  480(1),  420(1). 


TABLE  2 


Errors  in  £P(E)  and  N(E)  for  Different 
Grain  Sizes  and  Rounding  Methods 


Vibrational  Frequency  Mode 1 s_ 


Frequencies  (cm- 

')  After  Rounding 

Exact 

Frequencies  (cm-1) 

G/cm_I 

= 200 

A 

3000(2) 

3000(2) 

3000(2) 

1500(3) 

1400(3) 

1600(2) ,1400(1) 

1000(2) 

1000(2) 

1000(2) 

900(2) 

800(2) 

1000(1), 800(1) 

500(2) 

400(2) 

600(1) ,400(1) 

400(1) 

400(1) 

400(1) 

Sums  of  States 


E/103  cm-1 

Exact  Sums 

Percentage 

Deviation^ 

10 

15 

20 

24.4 

A 

5.394  x 105 
1.286  x 106 
1.529  x 10B 
9.363  x 10° 

Densities  of  States 

B 

99.9 

108.9 

114.9 
118.3 

C 

5.8 

4.4 

3.5 

2.9 

E/103  cm-1 

Exact  Densities 

Percentage 

__  . , d 

Deviation 

A 

B 

c 

10 

372.0 

102.2 

5.0 

15 

6.808  x 103 

109.9 

3.3 

20 

6.603  x 101- 

114.9 

2.5 

24.4 

3.523  x 10’ 

115.6 

0.5 

nSame  as  given  by  HH. 

^All  frequencies  rounded  downward  or  unchanged, 
c 

Rounding  scheme  as  given  in  text. 

^It  is  important  to  realize  that  these  "errors"  are  entirely  due  to  the  use  of 
frequencies  rounded  to  multiples  of  a grain  size  that  is  much  larger  than  as 
would  be  used  for  accurate  computations. 


Percentage  Deviations  in  Sums  of  States  for  Three 
Anharmonic  Oscillator  Models  as  a Function  of  Grain  Size 


E/103  cm-1 

£P(E) 

Percentage  Deviation  from  G/cm 

1 = 10  Model 

G/cm-1  =10 

G/cm-1  = 20 

G/cm-1  = 50 

G/cm-1  = 100 

Model  F 

5 

6155 

0.47 

1.8 

7.05 

10 

6.032  x 103 

0.66 

1.92 

5.95 

15 

1.556  x 107 

0.45 

1.29 

4.05 

20 

1.985  x 10s 

0.35 

1.06 

3.07 

24.4 

1.291  x 109 

0.23 

0.93 

2.48 

Model  G 

5 

6068 

0.00 

0.00 

0.13 

10 

5.484  x Xv. 

0.00 

0.29 

1.26 

15 

1.331  x 107 

0.00 

0.60 

1.65 

20 

1.613  x 10s 

0.00 

0.62 

1.61 

24.4 

1.003  x 109 

0.00 

0.60 

1.60 

Model  H 

5 ' 

6068 

0.00 

0.00 

1.45 

10 

5.503  x 10s 

0.09 

0.71 

3.71 

15 

1.349  x 107 

0.15 

0.89 

3.41 

20 

1.651  x 10s 

0.12 

0.85 

2.79 

24.4 

1.036  x 109 

0.10 

0.68 

2.32 

Models  F,  G,  and  H correspond  to  models  given  in  Reference  [4]  as  7,  9,  and 
10,  respectively.  These  models  use  the  same  six  fundamental  frequencies  given 
in  Table  1,  with  diagonal  anharmonicity  constants  (cm-1)  of 
Model  F * (-80,  -12,  -8,  -7,  0,  0) 

Model  G = (-50,  -5,  0,  0,  0,  0) 

Model  H = (-80,  -5,  0,  0,  0,  0). 
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FIGURES  la  and  lb 

Plots  of  P (E)  Versus  E for  120  Frequency  Harmonic  Oscillator  Model 

Fundamental  frequencies  in  cm-1  (degeneracies)  = 3000(20),  1500(30),  1000(20), 
900(20),  500(20),  400(10)  with  G = 100  cm-1.  Since  all  frequencies  are  exact 
multiples  of  G,  results  are  exact.  Points  denote  P(E)  values,  lines  are  drawn 
for  clarity  only. 


